
#include <math.h>

#include <util/common_macros.h>
#include <util/math_macros.h>
#include <util/msg.h>

#include "../../matrix/matrix.h"
#include "../image.h"
#include "../transforms/integral.h"
#include "threshold.h"

/**
 * \brief Apply a threshold to each element in array.
 * 
 * This function apply a threshold to each element in the input array.<br>
 * <br>
 * 
 * @param Img input array.
 * @param th reference threshold.
 * @param minVal botton value apply to elements less or equal to \c th.
 * @param maxVal top value apply to elements grater than \c th.
 * @param type define how the thresholding will be performed.
 * @param BW the threshold matrix result. If this parameter is NULL, the threshold is apply to input array \c Img.
 * 
 * \note
 * \arg If \c *BW is NULL, \c *BW is dynamically allocated in this function.<br>
 * \arg If \c BW is NULL, the threshold is applied at \c Img.<br>
 * 
 * \sa ImpThresholdType
 */
void imp_image_threshold(ImpMatrix * Img, real th, real minVal, real maxVal, ImpThresholdType type, ImpMatrix ** BW) {

	long int i;
	real * img;
	real * thresh;
	long int size;

	_exit_if_fail(Img);

	size = Img->size;
	img = Img->data;

	if (BW != NULL) {
		if (*BW == NULL) (*BW) = imp_matrix_new(Img->rows, Img->columns);
		thresh = (*BW)->data;
	} else {
		thresh = Img->data;
	}

	switch (type) {
		case IMP_THRESHOLD_BINARY:
			for (i = 0; i < size; i++) thresh[i] = (img[i] > th) ? maxVal : minVal;
			break;

		case IMP_THRESHOLD_BINARY_INV:
			for (i = 0; i < size; i++) thresh[i] = (img[i] > th) ? minVal : maxVal;
			break;

		case IMP_THRESHOLD_TRUNCATE_UP:
			for (i = 0; i < size; i++) thresh[i] = (img[i] > th) ? maxVal : img[i];
			break;

		case IMP_THRESHOLD_TRUNCATE_DOWN:
			for (i = 0; i < size; i++) thresh[i] = (img[i] > th) ? img[i] : minVal;
			break;

		case IMP_THRESHOLD_TO_ZERO:
			for (i = 0; i < size; i++) thresh[i] = (img[i] > th) ? img[i] : 0;
			break;

		case IMP_THRESHOLD_TO_ZERO_INV:
			for (i = 0; i < size; i++) thresh[i] = (img[i] > th) ? 0 : img[i];
			break;

		default:
			msg_error("Invalid threshold type!\n");
	}
}

/**
 * \brief This function calculates a binary image using Adaptive Thresholding.
 * 
 * This function support the methods defined by ImpThresholdAdaptiveMethod.
 * 
 * \see ImpThresholdAdaptiveMethod
 *  
 * @param Img input image.
 * @param type adaptive thresholding technique.
 * @param BW output image.
 */
void imp_image_threshold_adaptive(ImpMatrix * Img, ImpThresholdAdaptiveMethod method, ImpMatrix ** BW) {
	switch (method) {

		case IMP_THRESHOLD_ADAPTIVE_SAUVOLA_FAST:
			imp_image_threshold_Sauvola_fast(Img, Img->columns/15, 0.2, BW);
			break;

		case IMP_THRESHOLD_ADAPTIVE_WELLNER_EXTENDED:
			imp_image_threshold_wellner_extended(Img, BW);
			break;

		default:
			msg_error("Invalid adaptive threshold type!\n");
	}
}

/**
 * \brief Extended Wellner's Adaptive Thresholding.
 * 
 * This function implements the Extended Wellner's Adaptive Thresholding.
 *  
 * \sa Bradley, G., Roth, G. (2007), "Adaptive Thresholding Using Integral Image", Journal of Graphics Tools, Vol.12, No.2, p.13-21.
 *  
 * @param Img
 * @param BW
 */
void imp_image_threshold_wellner_extended(ImpMatrix * Img, ImpMatrix ** BW) {

	ImpMatrix * IntImg = NULL;
	real * img;
	real * imgi;
	real * thresh;
	real area;

	//int s = Img->columns / 8;
	int s = Img->columns / 6;
	int s2 = s / 2;
	int shift = s2 * Img->columns;
	int count = (s2 + 1) * (s2 + 1) * 4;
	real t = 0.15;

	img = Img->data;

	if (*BW == NULL) (*BW) = imp_matrix_new(Img->rows, Img->columns);
	thresh = (*BW)->data;

	imp_image_integral(Img, &IntImg);
	imgi = IntImg->data;

	int r, c, _r, i;
	int rows, columns;
	int i1, i2, i3, i4;

	rows = Img->rows;
	columns = Img->columns;

	for (r = s2, _r = r * columns; r < rows - s2; r++, _r += columns) {
		for (c = s2; c < columns - s2; c++) {
			i = _r + c;

			i1 = i - shift - s2;
			i2 = i - shift + s2;
			i3 = i + shift + s2;
			i4 = i + shift - s2;

			area = imgi[i3] - imgi[i2] - imgi[i4] + imgi[i1];
			if (img[i] * count <= area * (1 - t)) {
				thresh[i] = 0.0;
			} else {
				thresh[i] = 255.0;
			}
		}
	}

	int x1, y1, x2, y2;
	int rows1, columns1;

	rows1 = rows - 1;
	columns1 = columns - 1;

	// <editor-fold defaultstate="collapsed" desc="Check Top Border">
	for (r = 0, _r = r * columns; r < s2; r++, _r += columns) {

		y1 = r - s2;
		y2 = r + s2;

		y1 = CLAMP(y1, 0, rows1);
		y2 = CLAMP(y2, 0, rows1);

		for (c = 0; c < columns; c++) {
			i = _r + c;

			x1 = c - s2;
			x2 = c + s2;

			x1 = CLAMP(x1, 0, columns1);
			x2 = CLAMP(x2, 0, columns1);

			count = (x2 - x1)*(y2 - y1);

			i1 = x1 + y1*columns;
			i2 = x2 + y1*columns;
			i3 = x2 + y2*columns;
			i4 = x1 + y2*columns;

			area = imgi[i3] - imgi[i2] - imgi[i4] + imgi[i1];

			if (img[i] * count <= area * (1 - t)) {
				thresh[i] = 0.0;
			} else {
				thresh[i] = 255.0;
			}

		}
	}
	// </editor-fold>

	// <editor-fold defaultstate="collapsed" desc="Check Bottom Border">
	for (r = rows - s2, _r = r * columns; r < rows; r++, _r += columns) {

		y1 = r - s2;
		y2 = r + s2;

		y1 = CLAMP(y1, 0, rows1);
		y2 = CLAMP(y2, 0, rows1);

		for (c = 0; c < columns; c++) {
			i = _r + c;

			x1 = c - s2;
			x2 = c + s2;

			x1 = CLAMP(x1, 0, columns1);
			x2 = CLAMP(x2, 0, columns1);

			count = (x2 - x1)*(y2 - y1);

			i1 = x1 + y1*columns;
			i2 = x2 + y1*columns;
			i3 = x2 + y2*columns;
			i4 = x1 + y2*columns;

			area = imgi[i3] - imgi[i2] - imgi[i4] + imgi[i1];

			if (img[i] * count <= area * (1 - t)) {
				thresh[i] = 0.0;
			} else {
				thresh[i] = 255.0;
			}

		}
	}// </editor-fold>

	// <editor-fold defaultstate="collapsed" desc="Check Left Border">
	for (r = s2, _r = r * columns; r < rows - s2; r++, _r += columns) {

		y1 = r - s2;
		y2 = r + s2;

		y1 = CLAMP(y1, 0, rows1);
		y2 = CLAMP(y2, 0, rows1);

		for (c = 0; c < s2; c++) {
			i = _r + c;

			x1 = c - s2;
			x2 = c + s2;

			x1 = CLAMP(x1, 0, columns1);
			x2 = CLAMP(x2, 0, columns1);

			count = (x2 - x1)*(y2 - y1);

			i1 = x1 + y1*columns;
			i2 = x2 + y1*columns;
			i3 = x2 + y2*columns;
			i4 = x1 + y2*columns;

			area = imgi[i3] - imgi[i2] - imgi[i4] + imgi[i1];

			if (img[i] * count <= area * (1 - t)) {
				thresh[i] = 0.0;
			} else {
				thresh[i] = 255.0;
			}

		}
	}// </editor-fold>

	// <editor-fold defaultstate="collapsed" desc="Check Right Border">
	for (r = s2, _r = r * columns; r < rows - s2; r++, _r += columns) {

		y1 = r - s2;
		y2 = r + s2;

		y1 = CLAMP(y1, 0, rows1);
		y2 = CLAMP(y2, 0, rows1);

		for (c = columns - s2; c < columns; c++) {
			i = _r + c;

			x1 = c - s2;
			x2 = c + s2;

			x1 = CLAMP(x1, 0, columns1);
			x2 = CLAMP(x2, 0, columns1);

			count = (x2 - x1)*(y2 - y1);

			i1 = x1 + y1*columns;
			i2 = x2 + y1*columns;
			i3 = x2 + y2*columns;
			i4 = x1 + y2*columns;

			area = imgi[i3] - imgi[i2] - imgi[i4] + imgi[i1];

			if (img[i] * count <= area * (1 - t)) {
				thresh[i] = 0.0;
			} else {
				thresh[i] = 255.0;
			}
		}
	}// </editor-fold>

	imp_object_delete(IntImg);
}

/**
 * \brief Fast Sauvola's Adptive Thresholding.
 * 
 * This function calculates a binary image using a Sauvola algorithm.
 * The Fast Adaptive Image Binarization uses Sauvola's test and Summed Area Tables.
 * 
 * \see Faisal Shafait, Daniel Keysers, Thomas Breuel, (2008), "Efficient Implementation of Local Adaptive Thresholding Techniques Using Integral Images", Proceedings of the 15th Document Recognition and Retrieval Conference, IS&T/SPIE International Symposium on Electronic Imaging, January 26-31, San Jose, CA, USA volume 6815.
 * \see Sauvola, J & Pietaksinen, M (2000), "Adaptive Document Image Binarization", Pattern Recognition 33(2): 225-236
 * 
 * @param Img
 * @param w
 * @param k
 * @param BW
 */
void imp_image_threshold_Sauvola_fast(ImpMatrix * Img, int w, real k, ImpMatrix ** BW) {

	ImpMatrix * IntImg = NULL;
	ImpMatrix * SqIntImg = NULL;
	real * img;
	real * imgi;
	real * imgi2;
	real * thresh;
	real area, area2;

	/* Adaptive threshold calculated by each window */
	real Th;
	/* Maximum value of standard deviation */
	real max_s = 128.0;

	/* mean */
	real u;
	/* variance */
	real s2;
	/* standard deviation */
	real s;

	/* number of neighbors */
	int nn = w / 2;
	w = (2 * nn) + 1;
	int shift = nn * Img->columns;
	/* window size */
	int ws = w*w;
	//ws = (2n+1)(2n+1) = 4n² + 4n + 1;

	img = Img->data;

	if (*BW == NULL) (*BW) = imp_matrix_new(Img->rows, Img->columns);
	thresh = (*BW)->data;

	imp_image_integral2(Img, &IntImg, &SqIntImg);
	imgi = IntImg->data;
	imgi2 = SqIntImg->data;

	int r, c, _r, i;
	int rows, columns;
	int i1, i2, i3, i4;

	rows = Img->rows;
	columns = Img->columns;

	for (r = nn, _r = r * columns; r < rows - nn; r++, _r += columns) {
		for (c = nn; c < columns - nn; c++) {
			i = _r + c;

			i1 = i - shift - nn;
			i2 = i - shift + nn;
			i3 = i + shift + nn;
			i4 = i + shift - nn;

			area = imgi[i3] - imgi[i2] - imgi[i4] + imgi[i1];
			area2 = imgi2[i3] - imgi2[i2] - imgi2[i4] + imgi2[i1];

			u = area / ws;
			s2 = (area2 - u * u) / ws;

			s = sqrt(s2);
			
			/* adaptive threshold */
			Th = u * (1 + (k * ((s / max_s) - 1)));

			if (img[i] > Th) {
				thresh[i] = 255.0;
			} else {
				thresh[i] = 0.0;

			}
		}
	}
	
	/* = Border treatment =================================================== */
	int x1, y1, x2, y2;
	int rows1, columns1;

	rows1 = rows - 1;
	columns1 = columns - 1;

	// <editor-fold defaultstate="collapsed" desc="Check Top Border">
	for (r = 0, _r = r * columns; r < nn; r++, _r += columns) {

		y1 = r - nn;
		y2 = r + nn;
		y1 = CLAMP(y1, 0, rows1);
		y2 = CLAMP(y2, 0, rows1);

		for (c = 0; c < columns; c++) {
			i = _r + c;

			x1 = c - nn;
			x2 = c + nn;
			x1 = CLAMP(x1, 0, columns1);
			x2 = CLAMP(x2, 0, columns1);

			ws = (x2 - x1)*(y2 - y1);

			i1 = x1 + y1*columns;
			i2 = x2 + y1*columns;
			i3 = x2 + y2*columns;
			i4 = x1 + y2*columns;

			area = imgi[i3] - imgi[i2] - imgi[i4] + imgi[i1];
			area2 = imgi2[i3] - imgi2[i2] - imgi2[i4] + imgi2[i1];

			u = area / ws;
			s2 = (area2 - u * u) / ws;

			s = sqrt(s2);
			
			/* adaptive threshold */
			Th = u * (1 + (k * ((s / max_s) - 1)));

			if (img[i] > Th) {
				thresh[i] = 255.0;
			} else {
				thresh[i] = 0.0;

			}
		}
	}
	// </editor-fold>

	// <editor-fold defaultstate="collapsed" desc="Check Bottom Border">
	for (r = rows - nn, _r = r * columns; r < rows; r++, _r += columns) {

		y1 = r - nn;
		y2 = r + nn;

		y1 = CLAMP(y1, 0, rows1);
		y2 = CLAMP(y2, 0, rows1);

		for (c = 0; c < columns; c++) {
			i = _r + c;

			x1 = c - nn;
			x2 = c + nn;

			x1 = CLAMP(x1, 0, columns1);
			x2 = CLAMP(x2, 0, columns1);

			ws = (x2 - x1)*(y2 - y1);

			i1 = x1 + y1*columns;
			i2 = x2 + y1*columns;
			i3 = x2 + y2*columns;
			i4 = x1 + y2*columns;

			area = imgi[i3] - imgi[i2] - imgi[i4] + imgi[i1];
			area2 = imgi2[i3] - imgi2[i2] - imgi2[i4] + imgi2[i1];

			u = area / ws;
			s2 = (area2 - u * u) / ws;

			s = sqrt(s2);
			
			/* adaptive threshold */
			Th = u * (1 + (k * ((s / max_s) - 1)));

			if (img[i] > Th) {
				thresh[i] = 255.0;
			} else {
				thresh[i] = 0.0;

			}
		}
	}// </editor-fold>

	// <editor-fold defaultstate="collapsed" desc="Check Left Border">
	for (r = nn, _r = r * columns; r < rows - nn; r++, _r += columns) {

		y1 = r - nn;
		y2 = r + nn;

		y1 = CLAMP(y1, 0, rows1);
		y2 = CLAMP(y2, 0, rows1);

		for (c = 0; c < nn; c++) {
			i = _r + c;

			x1 = c - nn;
			x2 = c + nn;

			x1 = CLAMP(x1, 0, columns1);
			x2 = CLAMP(x2, 0, columns1);

			ws = (x2 - x1)*(y2 - y1);

			i1 = x1 + y1*columns;
			i2 = x2 + y1*columns;
			i3 = x2 + y2*columns;
			i4 = x1 + y2*columns;

			area = imgi[i3] - imgi[i2] - imgi[i4] + imgi[i1];
			area2 = imgi2[i3] - imgi2[i2] - imgi2[i4] + imgi2[i1];

			u = area / ws;
			s2 = (area2 - u * u) / ws;

			s = sqrt(s2);
			
			/* adaptive threshold */
			Th = u * (1 + (k * ((s / max_s) - 1)));

			if (img[i] > Th) {
				thresh[i] = 255.0;
			} else {
				thresh[i] = 0.0;

			}

		}
	}// </editor-fold>

	// <editor-fold defaultstate="collapsed" desc="Check Right Border">
	for (r = nn, _r = r * columns; r < rows - nn; r++, _r += columns) {

		y1 = r - nn;
		y2 = r + nn;

		y1 = CLAMP(y1, 0, rows1);
		y2 = CLAMP(y2, 0, rows1);

		for (c = columns - nn; c < columns; c++) {
			i = _r + c;

			x1 = c - nn;
			x2 = c + nn;

			x1 = CLAMP(x1, 0, columns1);
			x2 = CLAMP(x2, 0, columns1);

			ws = (x2 - x1)*(y2 - y1);

			i1 = x1 + y1*columns;
			i2 = x2 + y1*columns;
			i3 = x2 + y2*columns;
			i4 = x1 + y2*columns;

			area = imgi[i3] - imgi[i2] - imgi[i4] + imgi[i1];
			area2 = imgi2[i3] - imgi2[i2] - imgi2[i4] + imgi2[i1];

			u = area / ws;
			s2 = (area2 - u * u) / ws;

			s = sqrt(s2);
			
			/* adaptive threshold */
			Th = u * (1 + (k * ((s / max_s) - 1)));

			if (img[i] > Th) {
				thresh[i] = 255.0;
			} else {
				thresh[i] = 0.0;

			}
		}
	}// </editor-fold>

	imp_object_delete(IntImg);
	imp_object_delete(SqIntImg);
}

void imp_matrix_threshold(ImpMatrix * Img, real th, real minVal, real maxVal, boolean inv) {
	long int i;
	real * img;
	long int size;

	_verify(Img);

	size = Img->size;
	img = Img->data;

	for (i = 0; i < size; i++) {
		if (inv) img[i] = (img[i] <= th) ? maxVal : minVal;
		else img[i] = (img[i] <= th) ? minVal : maxVal;
	}
}

/*
void imp_image_threshold(ImpImage * Gray, unsigned char th, unsigned char minVal, unsigned char maxVal, int inv) {
	long int i;
	unsigned char * gray;
	long int size;

	_verify(Gray);

	size = Gray->size;
	gray = Gray->data;

	for (i = 0; i < size; i++) {
		if (inv) gray[i] = (gray[i] <= th) ? maxVal : minVal;
		else gray[i] = (gray[i] <= th) ? minVal : maxVal;
	}
}
 */


